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ABSTRACT 

We report Warm Spitzer full-orbit phase observations of WASP-12b at 3.6 and 4.5 //m. This 
extremely inflated hot Jupiter is thought to be overflowing its Roche lobe, undergoing mass loss, 
accretion onto its host star, and has been claimed to have a C/0 ratio in excess of unity. We are 
able to measure the transit depths, eclipse depths, thermal and ellipsoidal phase variations at both 
wavelengths. The large amplitude phase variations, combined with the planet's previously-measured 
day-side spectral energy distribution, is indicative of non-zero Bond albedo and very poor day-night 
heat redistribution. The transit depths in the mid-infrared — (i?p/i?*)^ = 0.0123(3) and 0.0111(3) at 
3.6 and 4.5 /im, respectively — indicate that the atmospheric opacity is greater at 3.6 than at 4.5 ^m, 
in disagreement with model predictions, irrespective of C/0 ratio. The secondary eclipse depths are 
consistent with previous studies: Fday/F^ = 0.0038(4) and 0.0039(3) at 3.6 and 4.5 iim, respectively. 
We do not detect ellipsoidal variations at 3.6 /im, but our parameter uncertainties — estimated via 
prayer-bead Monte Carlo — keep this non-detection consistent with model predictions. At 4.5 /im, on 
the other hand, we detect ellipsoidal variations that are much stronger than predicted. If interpreted 
as a geometric effect due to the planet's elongated shape, these variations imply a 3:2 ratio for the 
planet's longest: shortest axes and a relatively bright day-night terminator. If we instead presume that 
the 4.5 /xm ellipsoidal variations are due to uncorrected systematic noise and we fix the amplitude of 
the variations to zero, the best fit 4.5 /im transit depth becomes commensurate with the 3.6 /im depth, 
within the uncertainties. The relative transit depths are then consistent with a Solar composition and 
short scale height at the terminator. Assuming zero ellipsoidal variations also yields a much deeper 
4.5 fim eclipse depth, consistent with a Solar composition and modest temperature inversion. We 
suggest future observations that could distinguish between these two scenarios. 



1. INTRODUCTION 

Thermal phase variations are a powerful way to con- 
strain the climate on exoplanets. Such observations 
have been made for non-transiting short-period planets 
(|Cowan et al.l[2007HCrossfield et al.l [2010). but are most 
potent when combined with transit and eclipse obser- 
vations for edge-on systems, because of the additional 
knowledge of the planet's i nclination, mass and radius 
(jKnutson et al.ll2007ll2009bl laf). Secondary eclipse depths 
provide a constraint on the planet's day-side tempera- 
ture. Thermal phase variations probe the day-night tem- 
perature contrast and hence the planet's heat redistribu- 
tion efficiency. If the observational cadence and signal- 
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to-noise ratio are sufficiently high, phase variations are 
also sensitive to the offset between the noon meridian and 
the planet's hottest local stellar time, hence constraining 
wind speed and direction. 

By considering eclipse depths at a variety of 
wavelengths for a sa mple of 24 transiting planets, 
iCowan &: Agoll ()2011bD estimated their day-side effective 
temperatures, hence placing a joint constraint on the 
Bond albedo and heat recirculation efficiency of these 
planets. That study found that typical hot Jupiters ex- 
hibit a variety of albedo/recirculation efficiencies, but 
planets with substellar equilibrium temperatures greater 
than To ~ 2700 K all seem to have lower albedo and/or 
recirculation efficiency. In other words, the hottest tran- 
siting giant planets have a qualitatively different climate 
than the merely hot Jupiters, but it is not known whether 
this is due to a difference in albedo, circulation, or both. 
Direct measurements of hot Jupiter geometric albedos 
from optical secondary eclipse observations span more 
than an order of magnitude and do not resolve this de- 
generacy. 

In this paper, we br eak the albedo-rec irculation degen- 
eracy for WASP-12b (|Hebb et al.ll2009f ). one of the very 
hottest known exoplanets, with a day-side temperature 
of ^ 3000 K: the amplitude of thermal phase variations 
is a direct measure of the planet's day-night tempera- 
ture contrast and hence heat transport efficiency. If the 
night-side temperature is high, then the planet's albedo 
must be exceedingly low to be consistent with its high 
day-side temperature. If, on the other hand, the night- 
side temperature is low, then the planet has an albedo 
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in the tens of percent. 

WASP- 12b has been a fascinating planet since its 
discovery. The discrepant timing of its secondary 
ecli pse indicated that the plan et had a slight eccentric- 
ity dL opcz-Moralcs ct al. 2010) , but subsequent eclipse 
(ICampo ct al. 2011; CroU ct al., i2011[) and radial veloc- 



ity (jHusnoo et al. 2011) observations have all but ruled 
this out. Nevertheless, the planet's short-period orbit 
(1.1 days; just outside its star's Roche limit) and in- 
flated radius (1.8 Rj) led to th e prediction tha t it is 
tidall y disto rted (Ragozzine fc W oU 2009; Lccon te et al.l 
l2011t lBuda|l2011[ ). and undergoing Roch e-lobe overflow 
followed by ac cretion onto its host star (|Li et al.l 120101 : 
iLai et al1l2010[ ). The put ative early ing r ess of an ultra- 
violet transit observed by iFossati et al.l (|20100 seems to 
support this prediction, but may also be explained in 
terms of a leading bow shock from material stream ing 
off the planet ( Vidotto et al.ll2010l: iLlama et al.ir2011[ ). 



More recentiv. iMadhusudhan et al.l (|2011[ ) used the 



wavelength - depen danc e of mid-infrared ecli pse depths of 
iCroU et al.l (|2011[ ) and lCampo et al.l (|201in to constrain 
the atmospheric composition of WASP-12b, and found it 
has a carbon to oxygen ratio (C/0) greater than unity, 
unlike Solar System planets, or the assumed composi- 
tion of extrasolar planets. Those findings rested heavily 
on the relative eclipse depths at 3.6 and 4.5 fim. Our 
observations of eclipses and transits at these wavebands 
should be able to reinforce or rule out the high C/0 sce- 
nario. 

2. OBSERVATIONS & REDUCTION 

We acquired ob servations of WA SP- 12 (spectral type 
F9V) with IRAC (iFazio et al.l l 2004D on the Spitzer Space 
Telescope (|Werner et al.ll2004l ) at 3.6 fim (2010 Nov 17- 
18) and 4.5 ^m (2010 Dec 11-12) as part of the Warm 
Mission. We used the sub-array mode and acquired im- 
ages every 2 s (1.92 s effective exposure time), observing 
the system for approximately 1.3 days at each waveband, 
from slightly before a secondary eclipse to shortly af- 
ter the following secondary eclipse. This yields 902 data 
cubes (64 frames of 32 x 32 pixels) at each waveband. 

Due to a scheduling error, we did not observe all of the 
second eclipse's egress at 3.6 ^m. This does not severely 
affect our science objectives since we simultaneously fit 
both eclipses at a given wavelength; even at 3.6 /.tm we 
have nearly two full eclipse lightcurves to work with. 

We use the BCD files and convert MJy/str to 
electron counts by multiplying the flux values by 
GAIN xEXPTIME / FLUXCONV, using parameter val- 
ues from the header of the fits files. We use the 
BMJD_OBS and FRAMTIME parameter values to com- 
pute the BJD time at the middle of every exposure. 

We start by considering the pixel-by-pixel time-series 
for each 64-frame data cube, replacing NaN's with the 
pixel's median over that data cube; if the entire time- 
series for a given pixel is bad, it is flagged as a bad pixel 
and ignored in the subsequent analysis. At 4.5 /^m, the 
first r ow of pixe ls (y = 0) is consistently bad. 

Do ming et al.l (|2011f) noted that Warm Spitzer sub- 
array data cubes exhibit a frame-dependent background 
flux. At 3.6 /im, the first and 58th frames are consistent 
background outliers (both high) , and there is a clear drop 
in background flux througho ut each data cube (see Fig- 
ure 5 of .Deming et al.|[2011|) . At 4.5 /Ltm, the same two 



frames (1 and 58) are outliers (high and low backgrounds, 
respectively) , and there is a slight increase in background 
flux throughout each data cube. We elect to ignore the 
I''' and 58**^ frames of each data cube (3% of our data). 
To correct for the gradual change in background flux, we 
perform initial background subtraction on each frame of 
the data cube using the IDL routine MMM and excluding 
the 16 central pixels of the detector (those closest to the 
star). 

We then perform a two-step sigma clipping on each 
pixel's time-series, replacing 4cr outliers by the pixel's 
median in that data cube. The sigma-clipping at the 
pixel level affects 0.028% and 0.035% of our science time- 
series data at 3.6 and 4.5 /xm, respectively. 

To determine the centroid of our target, we flrst iden- 
tify the brightest of the central 16 pixels in each frame, 
then fit a 2D Gaussian to the 7x7 pixel box centered on 
this brightest pixel using the IDL routine GAUSS2DFItE3 

We perform aperture photometry on the individual 
frames of the data cubes using the IDL routine APER 
with a sky annulu s of 7-12 pixels in radius (as did 
ICampo et al.l[20Tll who observed the same system with 
the same instrument). Since we perform an initial back- 
ground subtraction early in our reduction pipeline (see 
above), our results are not very sensitive to changes in the 
sky annulus. We verified that nudging the inner/outer 
radius of the annulus by 1-2 pixels does not significantly 
affect the goodness-of-fit or astrophysical parameters. 

Bad pixels are ignored in the background estimation; 
images with a bad pixel within the aperture are ignored. 
To determine the optimal aperture for our analysis, we 
re-run our entire data reduction and analysis pipeline 
for a range of apertures, from 1.5 to 5.0 pixels, in in- 
crements of 0.5 pixel. We find that for both wavebands, 
the root-mean-squared scatter in the residuals are min- 
imized for an aperture of 2.5 pixels, which we therefore 
adopt for the remainder of our study. This is smaller 
than the apertures of 3.75 and 4.0 pixels (for 3.6 and 4 5 
fim images, respectively) used bv ICampo et al.l ()201lD . 
While using a larger aperture might reduce the photon- 
counting uncertainty, a small aperture makes it easier 
to correct for our dominant source of systematic uncer- 
tainty, the intra-pixel sensitivity variations described in 
Section 3.2. Since the Spitzer heater cycling was different 
for the two observing campaigns (see Section 3.2), it is 
possibl e that the nature of this systematic was different 
for the ICampo et al.l (|2011t ) observations. 

Lastly, we perform an iterative Aa clipping on the flux 
time-series, removing any outliers. This affects 0.01% 
and 0.02% of the science time-series data at 3.6 and 4.5 
fim^ respectively. The 3.6 /im cut is more generous be- 
cause of the greater systematic flux variations, as de- 
scribed below. 

The raw photometry is shown in Figure |TJ (Note that 
we perform all of our analysis on the unbinned data, but 
we bin the data for plotting.) The transits are easy to 
spot in the middle of the observations. The eclipses and 
phase variations are also visible by eye at 4.5 /xm. For 
the 3.6 fj,m light curve, the flrst eclipse and the phase 



We follow lAgol et aTl 120101 ). who compared many ccntroid- 
ing algorithms and found this one to be optimal. Using flux- 
weighted centroiding instead of PSF-fitting results in slightly worse 
X^, commensurate correlated noise as measured using /3 (see first 
Section 4.1), and consistent astrophysical parameters. 
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variations are difficult to distinguish from detector sys- 
tematics by eye. We estimate the system flux in mJy by 
converting back to MJy/str and using the pixel scale pa- 
rameter values, PXSCALl and PXSCAL2. Our system 
fluxes — 23.0(5) mJy and 14.7(f) mJy, at 3.6 and 4.5 /im, 
re spectively — are app roximately 10% lower than those 
of lCampo etaLl (12011]) . even when we adopt their larger 
apertures. iFazio et al.l ()2004D expected the absolute pho- 
tometric precision of IRAC to be better than 10%. 

As part of our Warm Spitzer observations, we also ob- 
tained mapping data in sub-array mode for WASP-12, at 
each 3.6 and 4.5 /im, immediately following the time se- 
ries described above. For the mapping data we acquired 
images every 0.4 s (0.36 s effective exposure time); we 
obtained 450 data cubes in each waveband. The purpose 
of these data was to map the central four pixels of the 
detector by scanning over them in 0.2-pixel and 0.1-pixel 
increments in the x- and t/-directions, respectively. 

The data reduction for the mapping data is identical to 
that for the science lightcurve, except that we stack the 
frames of each data cube using a pixel-by-pixel median. 
Aperture photometry is performed on these 450 stacked 
images rather than on the individual frames. This is nec- 
essary because of the shorter integration times for these 
data. (We also tried this data reduction — performing 
photometry cube-by-cube rather than frame-by-frame — 
on our phase curve data. Our best-fit models parame- 
ters were consistent using this approach, but the detector 
systematics were harder to correct for, leading to larger 
parameter uncertainties.) 

The pixel-level sigma-clipping affects 0.01% of both the 
3.6 and 4.5 /xm mapping data. There are no outliers in 
the flux time series for the mapping data. 

3. MODEL 

Our model has 9 free astrophysical parameters, plus 
up to 11 free parameters to characterize the detector re- 
sponse. The model parameters are listed in Table [1] and 
described below. 

3.1. Astrophysical Model 

Occultations: Transits are modeled u sing the IDL im- 
plementation of iMandel fc Agoll ()2002|) with fixed non- 
linear limb-darkening. To determin e the limb-da. rkening 
coefficients, we fit a four-parameter IClaretl ()2000[) model 
to a Kurucz stellar model wit h \Fe/H] = 0.3, Toff = 
6250 K, and log^ = 4.5 (jKumci IHll, 120(151 ). We 
set the eccentricit y to zero, and fix t he orb ital period 
to the value from iMacieiewski et all (|2011f) . The im- 
pact parameter, b, and the geometrical factor, a/R^,, 
are allowed to vary freely. Eclipses are modeled using 
the uniform-disk version of the Mandel & Agol (2002) 
expressions, and both eclipses are modeled simultane- 
ously with the same parameters. This means that we 
do not look for eclipse depth variability (searches for 
suc h variability have so far only resulted in upper lim- 
its: lAgol et alfeoiOl: iKnutson et all 1201 1[ ). We account 
for light travel time within the system, but this is only 
a matter of 22.5 s and does not affect our analysis. By 
the same to ken, we neglect eclip s e mapping effects for 
the planet (I Williams et al.l 120061: iRauscher et"all 120071: 
lAgol et al.ll2010[ ). since we are insensitive to the result- 
ing offset in eclipse time of less than a minute, let alone 
the detailed ingress/egress morphology. 



TABLE 1 
Model Variables 



Name Symbol 



Stellar Flux F* 

Orbital Period" P 

Impact Parameter b 

Geometric Factor a/i?. 

Time of Transit (BMJD) to 
Area Ratio (Rp/Rt)'^ 

Mean Planet Flux [Fp/F*) 

Thermal Phase Amplitude A^\\crui 

Phase Offset Omax 

Ellipsoidal Amplitude ^ellips 

t linear'' mt 

X Linear'^ ai 

y Linear"^ 6i 

X Quadratic'^ 02 

y Quadratic"^ b% 

X Cubic'^ 03 

y Cubic'^ 63 

X Quartic'^ 04 

y Quartic'^ 64 

X Quintic'^ as 

y Quintic'^ 65 



^ We fix the orbital period to the value from 
IMacieiewski et al.l I I2011I) . 

^ This parameter is only used when fit- 
ting occultations independently of the rest 
of the lightcurves. 

^ These parameters are only used in the 
polynomial IPSV fits described in Sec- 
tion [3TT3] 

Diurnal Phases: The planet's temperature and hence 
brightness vary as a function of local stellar time. This 
inhomogenous intensity is modeled with three parame- 
ters: the orbit-averaged planet/star flux ratio, (Fp/F^,), 
the semi-amplitude of thermal phase variations, Athcrm, 
and the offset of the phase peak from superior con- 
junction, Ofmax (ttmax < corrcsponds to a peak prior 
to superior conjunction and therefore to an eastward- 
advected hotspot and super-rotating winds). The phase 
variations have a sinusoidal shape, corresponding to a 
sinusoidal longi tudinal brightness profile for the planet 
(jCowan fc Agoli,2008, ')r'l Note that in the limit of poor 
recirculation, the longitudinal temperature profile should 
be more akin to a half-sine (uniformly dark on the night- 
side), leading to thermal phase variations similar to the 
Lambert phase function: broader minimum, briefer max- 
imum. We neglect reflected star light, which does not 
contribute appreciably at these wavelengths. 

Ellipsoidal Variations: Because of the planet's small 
semi-major axis and inflated radius, it is likely that it — 
and possibly its host star — are ellipsoidal in shape rather 
than spherical. This leads to changes in cross-sectional 
area throughout the planet's orbit. To good approxima- 
tion, these variations are sinusoidal with a period half 
of the orbital period; the maxima occur at quadrature, 
when we are seeing the long axes of the star and planet, 
and minima at superior and inferior conjunction, when 

In general, the offset between thermal phase maximum and 
superior conjunction is not the same as the offset of the hottest 
longitude of the planet with respect to the sub-stellar meridian. 
For realistic longitudinal temperature profiles, the observed offset 
in the lightcurve is si gnificantly greater than the hot spot offset 
llCowan Ik AgoM2011ai ). In the case of a first-order sinusoidal phase 
curve, however, the two offsets are one and the same. 
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where Athcrm and ^oiiips are the semi-amplitudes of diur- 
nal and ellipsoidal phase variations, respectively, and a 
is the phase angle (a = at superior conjunction, a — n 
at inferior conjunction). 

The secondary eclipse depth is related to the model 
variables by: 



= (^/ + ^thcrm COSftniax " 



ipsi 



Fig. 1. — Raw photometry in electron counts at 3.6 /^m (top) and 4.5 /im (bottom). The data have been binned for ease of viewing 

w e are see in g the sh ort axes of the t wo bo diesF^ 

iLi et all (poTot ). iLeconte et all (|20lTI ) and iBudail 
([2011') aU predict that the projected area of WASP-12b 
should vary by approximately 10% (peak-to-trough) due 
to its prolate shape, whether it is modeled as a pro- 
late ellipsoid or a partially-filled Roche Lobe. Given the 
mid-infrared planet/ star flux ratio of Fp/F^, « 4 x 10"'^ 
(jCampo et al.|[20Tll) . we expect to see ellipsoidal varia- 
tions in the planet at the level of AF/F^, w 4 x 10^'* (or 
a semi-amplitude of 2 x 10"^). 

The presence of a massive companion should also pro- 
duce ellipsoidal variations in the sta r, as seen at optica l 
wavelengths in the system HA T-P-7 ( Welsh et al."201()l). 
Using the expressions given in iFaigler fc Mazeh ( ,201l| , 
we estimate the semi-amplitude of these variations to be 
~ 4 X 10~^ at all wavelengths. We therefore expect that 
at thermal wavelengths, the ellipsoidal variations of the 
system should be dominated by the shape of the planet 
and not that of its host star. 

We experimented with different functional forms for 
the planet's phase variations in an effort to reduce corre- 
lations between astrophysical variables. Our best model 
in this regard is: 



(2) 



and to first order the associated uncertainties can be pro- 
pogated as: 



maxO-Athern 



^thcrm 



Sin^ Ofn 



(3) 



^ = (^)-|-AthormCOs(Q!-a„ 



-AompsCos(2Q;), (1) 



This is a common ap proximation for ellipsoidal variations 
(e.g.. FF^igier fc Mazei]ir2011l ). but we discuss the exact expression 
in Section 5.3.1. 



In practice, this is an overestimate of uncertainty, be- 
cause (Fp/F*) and Athcrm are anti-correlated. 

3.2. Detector Response Model 

IRAC channels 1 and 2 exhibit well-known intra- 
pixel sensitivity variations (IPSVs): photons hitting 
certain parts of a pixel lead to more electron counts 
than photons hitting other parts of the same pixel 
(e.g., iCharbonneau et aU l2005t iMorales-Calderon et al.l 
l2006f ) F^ In general, the sensitivity to photons is lowest 
near the edges of a pixel and greatest near its center (see 

This is entirely different than mier-pixel sensitivity variations, 
which should have been largely calibrated out by fiat-fielding. In 
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bottom panels of Figure [5]). On its own, IPSVs would 
not be a problem for precision time-resolved relative pho- 
tometry. But over the course of observations, the point 
spread function (PSF) of the target star moves on the 
detector. Even though the PSF spans many pixels, the 
IPSVs do not average out, because most of the flux falls 
in the PSF core. 

Since they are ultimately caused by changes in the PSF 
position, attempts to correct for IPSVs rely on accurate 
centroiding (described in Section 2). The centroiding is 
shown in Figure [5] for 3.6 fim (left) and 4.5 /im (right). 
The top panels show the centroid jitter and drift over 
the course of the observations; the second panels show 
the two-dimensional path of the centroid; the bottom 
panels show the intra-pixel sensitivity map of the four 
cen tral pixels of t he de tector, constructed by applying 
the Ballar d et al.l ()2010( ) point-by-point decorrelation to 
our mapping data (see Section r3.2.2p . 

Both X and y centroids exhibit fast jitter (period of 
roughly half an hour) with peak-trough amplitude of 
0.05-0.1 pixels. This is the same jitter that used to have 
a period of roughly an hour: it is related to the heater 
cycling on Spitzer. The cycling frequency was doubled in 
fall 2010, which doubled the centroid jitter frequency and 
roughly halved the amplitude of the jitter 13 The smaller 
amplitude of the jitter is undoubtedly an improvement, 
while the higher frequency may or may not be a nui- 
sance, depending on the duration of ingress/egress for a 
given planet. In our data, the 3.6 /im flux exhibits clear 
1% peak-to-trough flux variations on the centroid-jitter 
timescale, while the 4.5 iim flux does not. 

There is also a longer-term centroid drift, which is 
greatest in the y-direction: 0.5 pixel of motion over the 
'-^l day observation at both 3.6 and 4.5 /zm. The x- 
direction shows almost no long-term drift (0.05 pixels, 
comparable to the faster jitter). Looking at the bot- 
tom panels of Figure [51 it is easy to understand why 
the IPSVs are worse at 3.6 /im than at 4.5 /im: at the 
shorter waveband, the PSF drifted up a steep slope in 
sensitivity from a pixel corner towards a center; at the 
longer waveband, the PSF contoured below a ridge in 
sensitivity. 

The telescope takes a few hours to settle after pointing 
at a new target, resulting in larger PSF excursions for 
the first few data cubes of each lightcurve. It is difficult 
to correct for IPSVs in poorly-sampled regions of the 
detector, so we remove the first 0.05 days of both the 
3.6 and 4.5 /tm lightcurves (3.8% of our data in each 
waveband). 

The crux of the data reduction process is correcting for 
IPSVs, because our astrophysical signals (eclipses and 
phase variations) have an expected amplitude of 0.4%, 
comparable to — or smaller than — these systematics. 
We tried a variety of techniques, of which we describe 
the most promising below. We first present two methods 
for removing the IPSVs before fitting our astrophysical 
model of the system. These techniques are attractive be- 
cause they allow one to produce a systematics-corrected 
lightcurve independent of any astrophysical model as- 
sumptions. We then present an attempt to simultane- 

any case, a scheme that corrects IPSVs impUcitly corrects inter- 
pixel sensitivity variations as well. 

ssc.spitzer.caltech.edu/warmmission/news/21oct2010memo.pdf 



ously fit the IPSV and the astrophysical brightness vari- 
ations. 

3.2.1. Gaussian Decorrelation 

We follow iBallard et all (|2010[ ) in using point-by-point 
positions and fluxes to generate an intra-pixel sensitiv- 
ity map with x and y Gaussian smoothing le ngths of 
and cTy, respectively. IStevenson et al.l (j2011|) recently 
introduced a similar — but faster — method using bilin- 
ear interpolation. Since we are not attempting to itera- 
tively fi t the astrophysical a nd IPSV model, we use the 
simpler IBallard et al 1 ()2010D method. We experimented 
with different smoothing lengths and chose the combina- 
tions that yield the smallest value for the final model 
fit; the astrophysical parameters are not very sensitive 
to changes in the smoothing lengt h. At 3.6 /im we us e 
(7^ = 0.017 and ay = 0.0043, as in IBallard et all (f2010l) : 
at 4.5 /(m we use ax = ay = 0.015. The resulting pixel 
maps are shown in the top panels of Figure [H 

We then divide the raw photometry by the weight func- 
tion and fit our astrophysical model. The second panels 
of Figure [3] show the corrected data with best-fit astro- 
physical model, and residuals. The bottom panels show 
the scatter in the residuals as a function of binning, along 
with a red line indicating the Gaussian-noise limit of 
root-mean-squared scatter scaling as Vn. The normal- 
ization of this theoretical curve is based on the Poisson 
error for our electron counts (see first the raw photom- 
etry in Figure [1} . The best- fit astrophysical parameters 
are listed in Tables [2] (3.6 /im) and [3] (4.5 /im). 

3.2.2. Mapping Data 

The purpose of the mapping data was to deliberately 
map the central four pixels of the detector by scanning 
over them in 0.2-pixel and 0.1-pixel increments in the x- 
and j/-directions, respectively. Since these observations 
are much shorter than the planet's orbital time, and were 
scheduled to avoid transits or eclipses, the changes in flux 
are in principle entirely due to the centroid position on 
the detector. 

In practice, the 3.6 /im lightcurve observations ended 
approximately 3.8 minutes before the end of eclipse 
egress, and the 3.6 /im mapping observations immedi- 
ately followed. We therefore remove the first 4 minutes 
(0.003 days) of the 3.6 /im mapping data. 

We use the mapping data centroids and fluxes to gen- 
erate a weight map at the locations of the science cen- 
troids, again following [Ballard et al.l (|20Tol ). We adopt 
larger smoothing kernels set by the Nyquist sampling fre- 
quency of the regularly-spaced mapping centroids (ax = 
0.1, fjj, — 0.05). We then use this weight function to 
correct the science lightcurve as above. 

Using the mapping data to generate pixel maps has 
the advantage that we are not self-calibrating our sci- 
ence data, and hence are not liable to throw the baby 
out with the bath water. On the other hand, it does 
not successfully remove the systematics: the model fits 
are far worse than either the Gaussian decorrelation dis- 
cussed above, or the polynomial models discussed below 
(the discrepancy in is a factor of ~ 10 at 3.6 /tm and 
~ 4 at 4.5 /im). This means that the IPSV must have fine 
spatial structure (as seen by IBallard eraI|[20Tol) . and/or 
some additional flux- or time-dependence. 
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Fig. 2. — Point-spread-function (PSF) centroid movement at 3.6 fim (left) and 4.5 fim (right). The top panels shows the jitter and drift of 
the PSF centroid. The second panels show the two-dimensional wander of the centroid. The bottom panels shows the intra-pixel sensitivity 
map for the centr al reg ions of the detector, constructed by applying the lBallard et aLl COlOl ) point-by-point decorrelation to our mapping 
data (see Section I3.2.2I I. with red arrows marking the approximate drift of the PSF over the course of our observations. The green lines 
show pixel edges. 





Fig. 3. — WASP- 12b at 3.6 /im (left) and 4.5 (right), where we have corrected for intra-pixel sensitivity variations (IPSV's) using the 
I Ballard et al.l 112010 ) point-by-point decorrelation, as described in Section 3.2.1. The upper panels show the IPSV map determined from our 
science data, with pixel edges shown in green. The second panels show the corrected data with best-fit astrophysical model, and residuals. 
The bottom panels show the scatter in the residuals as a function of binning; the red line shows the photon noise limit; the vertical dotted 
line denotes the timescale of ingress/egress. 



3.2.3. Polynomial IPSV Model 

Here we model the intra-pixel sensitivity variations as a 
polynomial in the centroid x and y. We simultaneously 
fit our astrophysical model and the IPSVs by treating 
the X and y centroids as independent variables in our 
function. We model the IPSV's as: 

J^obs ^ 1 + Er=i [a^ix-xy + b,{y-yY] 



where i^obs is the observed flux, Fastro is the astrophysi- 
cal model, and x and y are the mean centroid positions. 
Formally, cross-terms are necessary to describe an arbi- 
trary two-dimensional function, but we find that includ- 
ing cross-terms does not significantly improve the or 
affect the astrophysical parameters. This is a testament 
to the fact that the bulk of the PSF motion is in the 
j/-direction. 

We experiment with polynomials up to sixth-order. To 
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test whether each additional pair of coefficients (one each 
for X and y) improved th e fit, we use t he Bayesian In- 
formation Criterion (BIC: ISchwarzlfTOTSl) . which imposes 
a penalty term on the for additional free parameters: 
BIC = + k\ivN , where k is the number of free pa- 
rameters and N k, 52000 is the number of data. Since 
ln(52000) sa 11, an additional parameter is acceptable if 
it improves the by at least 11. The BIC for our 3.6 
/xm data improved with the addition of parameters up to 
and including fifth order. The BIC for our 4.5 /im data 
was not improved by the addition of terms beyond cubic. 

Unlike some previous studies, we do not include a lin- 
ear ramp in time. Since the bulk of the PSF motion is a 
monotonic drift in the y-direction, we found the ramp in 
time to be highly correlated with the linear and quadratic 
coefficients of the y sensitivity, and the fit was not sig- 
nificantly improved 

We show the resulting fits in Figures |4] and El For each 
figure, the top panel shows: the systematics-corrected 
lightcurve and best-fit astrophysical model (top inset), 
the residuals after subtraction of the best-fit thermal 
phases, along with the best-fit ellipsoidal variations 
model (middle inset), and the residuals after removing 
ellipsoidal variations (bottom inset). (Ellipsoidal vari- 
ations of the planet do not affect in-eclipse data since 
the planet is hidden from view; hence we remove the 
in-eclipse data from this panel.) The bottom- left panel 
shows the weight function used to correct the data. The 
bottom-right panel shows the scatter in the residuals as a 
function of binning; the red line shows the photon noise 
limit; the vertical dotted line denotes the timescale of 
ingress/egress. 

4. MODEL FITTING & ERROR ANALYSIS 

We use the IDL implementation of a Levenberg- 
Marquardt (L-M) gradient descent routine, MPFITFUN, 
to find the best-fit model parameters. The covariance 
matrix of the model parameters provides a first guess at 
the parameter uncertainties. 

L-M or Markov Chain Monte Carlo (MCMC) error es- 
timation depend on the photometric uncertainties: larger 
error bars on the data lead to larger uncertainties on 
the model parameters. For the initial fits, we optimisti- 
cally set the error bars on our data at the Poisson limit, 
1 / V -^counts; this means that the reduced x^ of our best- 
fit model fits is somewhat larger than 1 (see Tables [5] 
and |3|), and it means that either L-M or MCMC will 
underestimate parameter uncertainties. To alleviate this 
problem, we then scale the data uncertainties to give a 
reduced x^j xIjj of unity (by multiplying uncertainties 
by the square-root of the best xlj)- In the present case, 
this entails inflating the error bars by a constant factor 
of 10-20%, depending on the waveband and model in 
question. Scaling photometric errors to obtain a reduced 
X^ of unity renders the x^ useless for comparing differ- 
ent models; we therefore quote the best reduced x^ prior 
to inflating the error bars. After adjusting the photo- 
metric uncertainties, we normalize the lightcurve to the 

We do include a linear ramp in time when performing isolated 
fits to transits or eclipses, since 1) those shorter time-series do not 
provide enough leverage to properly fit the IPSV, and 2) the linear 
ramp can act as a proxy for the thermal phase variations, which 
are not explicitly fit for in these cases. 



in-eclipse (star-only) value, and fix to unity in order 
to avoid correlations between and (Fp/F^) in the final 
fits. 

It is well known that simply using the covariance ma- 
trix from the L-M fit does not provide a robust error 
estimate, so we estimate parameter uncertainties using 
a variety of other techniques. We experimented with 
MCMC and Bootstrap Monte Carlo error estimates and 
found them to be slightly larger than — but compara- 
ble to — those from the L-M. We found that our most 
conservative error estimates (typically larger by a factor 
of 2 or more) are obtained by considering the residuals 
of our best-fit model: either binning of residuals, or re- 
sampling of residuals using a prayer-bead Monte Carlo. 
Throughout the manuscript we always adopt the largest 
uncertainty for a given parameter, but it is worth noting 
that even our most conservative error estimates may still 
be 15-30% small er than what one wou ld obtain with a 
wavelet analysis ([Carter fc Winnll2009[ ). 

4.1. Binning of Residuals 

iPont et all ((gbOG) proposed a simple method to ac- 
count for red noise by considering how the scatter in 
residuals decreases with bin size (bottom panels of Fig- 
ures 3, 4 and 5). At the left end of the plots, where we are 
considering point-to-point scatter in the residuals, the 
scatter is only 10-20% greater than the Poisson count- 
ing limit shown in red (cx ^/ M/N(M — 1) « \^N, where 
N is the number of observations per bin, and M is the 
number of bins). But the observed scatter does not fol- 
low the theoretical relation as the data are binned. The 
most important timescale for transit and eclipse param- 
eter estimation is the duration of ingress/egress, which 
we denote by a vertical dotted line in those panels (21 
minutes for WASP- 12b). The scatter on this timescale 
determines the accuracy we can expect to achieve for 
transit or e clipse depths . 

Following IWinn et al.l (|2007l . I2008| ) . we define the fac- 
tor /3 as the actual scatter (black line) divided by the 
theoretical Poisson limit (red line) on the 21 minute 
timescale (vertical dotted black line); our residuals have 
/3 — 2-3. To account for red noise in parameter uncer- 
tainties, we simply infiate the L-M parameter uncertain- 
ties by /3. (Inflating the photometric error bars by /3 
and recomputing the covariance matrix using L-M takes 
longer and produces slightly smaller parameter uncer- 
tainties.) The binning of residuals method turns out to 
be the most conservative error estimate for transit- and 
eclipse-specific model parameters (6, a/i?», (Rp/R^)^, 
etc.). 

4.2. Boot-Strap & Prayer-Bead Monte Carlo 

We also estimate parameter uncertainties using two re- 
sampling techniques: boot-strap and prayer-bead Monte 
Carlos. Both of these techniques use the scatter in the 
residuals of our best fit model as an estimate of photo- 
metric uncertainty. In both cases the residuals are shifted 
in time and added back to the best fit model to produce a 
new instance of the lightcurve, which is then fit using the 
L-MI3 The standard deviation in the sequence of model 

-^^ Since we use L-M to find best-fit solutions, one might think 
that the prayer-bead and boot-strap analyses implicitly depend on 



Thermal Phase Variations of WASP-12b 



9 



TABLE 2 

3.6 MICRON Parameters 



Calibration Method Xr 


b 


a/R, 




Fday/F* 


2^thcrm 




^cllips 


Point-by-Point Decorrelation 1.392" 
''Polynomial in x and y 1.384° 


0.5(1) 
0.3(2) 


2.8(2) 
3.1(2) 


0.0125(4) 
0.0123(3) 


0.0038(4) 
0.0033(4) 


0.0004(3) 
0.0038(6) 


0(29)° 
-53(7)° 


1(1) X 10-4 
1(1) X 10-4 



^ When comparing these values, it is worth remembering that with approximately 52,000 degrees of freedom, a model is 
signifieantly better if it improves x%, ^t least 0.004. 



Our fidueial analysis. 



TABLE 3 
4.5 MICRON Parameters 



Calibration Method Xr b a/R, {Rp/R,)'^ Fday/F* 2Athorm ^max -4enips 

Point-by-Point Decorrelation 1.326" 0.5(1) 2.9(2) 0.0112(4) 0.0039(3) 0.0019(3) -12(6)° 1.1(1) X IQ-^ 
''Polynomial in x and y 1.324° 0.5(1) 2.9(2) 0.0111(4) 0.0039(3) 0.0040(3) -16(4)° 1.2(2) X IQ-^ 

^ When eomparing these values, it is worth remembering that with approximately 52,000 degrees of freedom, a model is 
signifieantly better if it reduces X/i i^y at least 0.004. 
^ Our fiducial analysis. 



parameters is our estimate of their Icr uncertainty. Note 
that — by construction — resamphng techniques cannot 
improve parameter estimates: the best-fit parameters are 
those determined by fitting the original time series. 

For the boot strap Monte Carlo, the residuals are 
randomly shuffled so — much like the L-M and MCMC 
techniques — the bootstrap error analysis is insensitive 
to the ordering of residuals. Such error estimates will 
therefore only be accurate insofar as residuals are uncor- 
related in time. Indeed, we found that error estimates 
from the bootstrap were comparable to those from the 
L-M or MCMC analyses. 

The prayer-bead analysis maintains the relative order- 
ing of the residuals and simply shifts them all by the same 
amount (wrapping around the start/end of the data), so 
that correlated noise present in the residuals is preserved. 
A prayer-bead analysis is not appropriate if the nature 
of the noise is expected to change throughout the obser- 
vations (e.g., the 8 /im ramp seen in cryogenic Spitzer 
observations). But there is no evidence for such changes 
in behavior in Warm Spitzer data in general, or in our 
time series in particular. The prayer-bead can only have 
as many iterations as there are data points, but this is not 
a problem for the current study given our approximately 
52,000 images; we run the prayer bead for 10,000 iter- 
ations with randomly chosen offsets and verify that the 
uncertainties have converged by comparing to 100- and 
1000-iteration prayer-bead MCs. It is also worth noting 
that with such a long data set, we are more likely to ob- 
serve rare instances of bad behavior in the detector, mak- 
ing the prayer-bead technique particularly conservative. 
Indeed, prayer-bead Monte Carlo provides the largest 
error bars for the phase variation parameters: ^therm, 

^niax5 ^ellips- 

5. RESULTS 

The two methods used to remove intra-pixel sensitiv- 
ity variations are fundamentally different. The Gaussian 

photometric error estimates. But in fact, the L-M algorithm settles 
on the same solution irrespective of error bars (within reason). The 
prayer bead and bootstrap parameter uncertainties are therefore 
effectively independent of the photometric uncertainties. 



decorrelation uses local information to correct the flux 
with no assumption about larger-scale trends; it is able 
to correct for small-scale variations in sensitivity, but re- 
quires very high densities of centroids. The polynomial 
flt assumes a functional form for the smoothly-varying 
sensitivity, but is better able to correct regions that are 
less-well sampled. It is not c lear which me t hod is better 
suited to a given data set. iBallard et all (|2010[) found 
the decorrelation to be better (in terms of y^) than the 
polynomial flt when analyzing their binned 4.5 /.tm time 
series. We run both methods on unbinned data and flnd 
that the two methods perform equally well at 4.5 /im, 
while at 3.6 /xm the polynomial fit is better. 

Given their very different underlying philosophies, it is 
encouraging that the two methods yield similar weight 
functions (compare the pixel maps in Figure 3 to those 
in Figures 4 and 5). In fact, the transit depths, eclipse 
depths, and ellipsoidal variations recovered by the two 
techniques are generally consistent. The thermal phase 
variation amplitude and offset differ significantly, how- 
ever. 

5.1. Transits 

Our values for the impact parameter and geometrical 
factor are broadly consistent with published values. Note 
that we simultaneously fit the transits and eclipses, and 
eclipses are notoriously bad at constraining these orbital 

parameters 

Combining our transit times for epochs 925 and 
947 (BMJD of 55,518.040 7(4) and 55 , 542.05 21(4), re- 
spectively) with those of | Chan et all (|2011| ) and the 
ephemeris of lMacieiewski et al.l ()2011t ). we obtain a BJD 
discovery epoch transit center of to = 2454508.9768(2) 
and period P — 1.0914207(4) days. Since our transits oc- 
curred slightly earlier than predicted, our best-fit orbital 
period is 3.5f7 shorter than that of iMaciejewski et ahl 
(|2011h . but it is notoriously difficult to compare tran- 
sit times at different wavelengths analyzed by different 

If we only consider the 0.3 days of data centered on each 
transit, we find impact parameters of 6 = 0.46 and 0.60 (larger 
than published values), and geometrical factors of a/ij, = 2.9 and 
2.7 (smaller than published values) at 3.6 and 4.5 /im, respectively. 
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Fig. 4. — WASP-12b at 3.6 /xm, where we have treated the IPSVs as a polynomial function in both x and y centroid, as described in 
Section 3.2.3. The top panel shows: the systematics-corrected lightcurve and best-fit astrophysical model (top inset), the residuals after 
subtracting the best- fit transit, eclipse and thermal phase model, along with the best-fit ellipsoidal variations model (middle inset), and 
the residuals after removing ellipsoidal variations (bottom inset). Ellipsoidal variations of the planet do not affect in-eclipse data since the 
planet is hidden from view; hence we remove the in-eclipse data from this panel. The bottom-left panel shows the weight function used to 
correct the data, with pixel edges shown in green. The bottom-right panel shows the scatter in the residuals as a function of binning; the 
red line shows the photon noise limit; the vertical dotted line denotes the timescale of ingress/egress. 
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groups because of subtle differences in star-spot coverag e, 
treatment of limb- darkening, etc. (jPesert et al.ll20lfbf) . 

More importantly, we are simultaneously fitting an en- 
tire orbit worth of data including many different astro- 
physical effects, while the optical transit data were fit 
on their own, irrespective of longer-term astrophysical 
trends. In order to make a more appropriate compar- 
ison, we try fitting the transits independently of the 
rest of the data (considering only those data within 0.15 
days of the transit center). This yields later transit 
centers , leading to an ep hemeris more consistent with 
that of iMacieiewski et all (poTTl ): to = 2454508.9767(2), 
P = 1.0914210(4). 

Because of WASP-12b's high temperature and inflated 
radius, the variations in transit depth with wavelength 
should be 3x larger than for "typical " hot Jupiters (using 
the scaling relation from Winn 2010). Ba sed on observa - 
tions and models of HD 189733b [Fortnev et al.l[2010l ). 
one might therefore expect the transit depth at 4.5 fj,m 
to be ~ 10"'^ deeper than at 3.6 /xm. This is borne 
out by the dot t ed black lin e in Figure 6, which shows a 
iBurrows et al.l (|2007l l2008| ) model transit spectrum as- 
suming Solar composition and a day-like temperature- 
pressure (T-P) profile. If the planet's terminator has a 
night-like T-P profile (shorter scale height; solid black 
line in Figure 0, the difference in transit depth could be 
as small as 2 x 10~^ — but always with a transit depth 
greater at 4.5 //m than at 3.6 /xm. 

In Figure 6 we compare our transit depths at 3.6 
and 4.5 /xm with previously p ublished o ptical values: 
0.0138(2) (B and z'-band , Hebb ^al.ll2009h . 0.01380(16) 
(R-band, 'Macieiewski e±aL| mi^, and 0.0125(4) (V- 
band, Chan et al. 2011), yielding a three-band transmis- 
sion spectrum of the planet. Curiously, the transit depth 
at 4.5 fiTa is con siderably shallower than a t 3.6 /im. If we 
adopt the larger lMacieiewski et al.l ()2011f ) optical transit 
depth, then our mid-IR transit depths could be indicative 
of hazes or optical absorbers in the atmospher e of WASP- 
12b, as has been inferred for HD 189733b (iPont et al.l 
[2008; Sing ct al. 2011). The larger radius at 3.6 ^m as 
compared to 4.5 fim indicates a higher atmospheric opac- 
ity at the shorter waveband, which is difficult to reconcile 
with current models. 

In an attempt to e xplain its peculiar eclipse spectrum, 
iDeming etHI (|2011h hypothesized that CoRoT-2b might 
have equal opacity at 3.6 and 4.5 fim (and lower opacity 
at 8 fXTo) due to a haze of — as yet unknown — /xm-sized 
particles. One may similarly explain the unusual WASP- 
12b transit spectrum in terms of a haze of slightly smaller 
particles, such that the opacity drops from 3.6 to 4.5 
/im. Given the large difference in transit depth between 
the two wavebands, this hypothesis also requires a large 
atmospheric scale height at the planet's terminator. 

5.2. Eclipses 

At 3.6 fiTa, the eclipse depth is ^ la lower using 
the polynomial fit as compared to the decorrelation. If 
we fit the two eclipses individually using the polyno- 
mial IPSV fit, we obtain depths of 0.0030 (highly corre- 
lated residuals) and 0.0038 (incomplete egress), respec- 
tively. Given the prior measurement of 0.0038(1) by 
iCampo et al.l ()2011l ). we adopt the larger value from the 
Gaussian decorrelation for our analysis (this choice does 
not significantly affect any of our conclusions) . 



At 4.5 ^m, we obtain comparable values and eclipse 
depths regardless of our treatment of sy stematics. In 
both cases they are consistent with the iCampo et al.l 
(|20ll value: 0.0038(2). 

In all cases our erro r estimates are somewhat larger 
than the Campo et al.l (]2011l ) estimates. We observed 
the same planetary system with the same instrument, 
so one expects the same eclipse depths and uncertain- 
ties. The only significant change in the detector is that 
our observations occurred after the fall 2010 change in 
heater cycling, as mentioned in Section 2. This means 
that the short-term telescope jitter for our observations 
has a period of 30 minutes rather than 1 hour, and the 
amplitude of the centroid excursions — and hence flux 
variations — is reduced by a factor of two (for reference, 
the ingress/egress time for WASP-12b is 21 minutes, and 
the total transit/eclipse duration is 3 hours). 

It is conceivable that the near coincidence between 
the centroid jitter half-period and the ingress/egress 
timescale leads to greater res i dual systematics in our 
data than in the ICampo et al.l ()2011f) time series. How- 
ever, our eclipse depths are based on two occultations, 
and we have a much longer baseline of observations to 
help us correct for detector systematics, characterize 
noise properties and estimate uncertainties (see Section 
4. 2). Our MCMC err or estimates are similar to those 
of lCampo et al.l (|2011f ). but residual-binning and prayer- 
bead analyses produce eclipse depth uncertainties more 
than 2x larger than the MCMC. This means that there 
is still red noise present in our residuals, and the prayer- 
bead analysis is a more realistic estimate of our param- 
eter uncertainties. 

5.2.1. Day-Side Emergent Spectrum 

If one assumes solar composition, the relative eclipse 
depths at 3.6 and 4.5 fiui can probe the temperature vs. 
pressure (T-P) profile of the planet. Water vapor ab- 
sorbs less at 3.6 than at 4.5 /im, so the shorter waveband 
will have a higher brightness temperature if the temper- 
ature is loca lly dropping with height, or vice versa (e.g., 
IBurrows fc"Ortonll2010f l. 

If the eclipse depth is measured at sufficiently many 
wavelengths, one may hope to simultaneously constrain 
a planet's atmos pheric composition and T-P profile (e.g., 
iMadhusudhMTfc Scager 20 0l). The high C/0 chem- 
istry invoked by M adhusudhan et al.l (|2011l ) was the re- 
sult of an abnormal ly low eclipse depth at 4.5 /im in the 
ICampo et al.l (j2011h data, which was interpreted as being 
due to CO absorption in the pl anet's relatively cool up - 
per atmosphere. More recently, iKopparapu et al.l (|201lD 
used a photochemical model to study the disequilibrium 
chemistry of WASP-12b, confirming that CO would be 
enhanced in a high C/0 composition atmosphere. 

In Figure [7l we compare the near- to mid-IR broad- 
band spectrum of WASP-12b to various 1- D radiative 
transfer models (|Burrows et all 120071 |2008[) . We vary 
the abundance of CO as a proxy for varying the C/0 
ratio. Our e clipse depths are consistent with those of 
ICampo et all ([201 1), so we still favor models with en- 
hanced CO (lOx Solar) a nd a weak inversion f or thi s 
planet, in agreement with iMadhusudhan et aTl ()2011l) . 
We also find that we can obtain an equally good fit to 
the data by reducing H2O to 1% Solar abundance and 
partitioning carbon evenly between CO and CH4. Given 
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Fig. 5. — WASP-12b at 4.5 /xm, where we have treated the IPSVs as a polynomial function in both x and y centroid, as described in 
Section 3.2.3. The top panel shows: the systematics-corrected lightcurve and best-fit astrophysical model (top inset), the residuals after 
subtracting the best- fit transit, eclipse and thermal phase model, along with the best-fit ellipsoidal variations model (middle inset), and 
the residuals after removing ellipsoidal variations (bottom inset). Ellipsoidal variations of the planet do not affect in-eclipse data since the 
planet is hidden from view; hence we remove the in-eclipse data from this panel. The bottom-left panel shows the weight function used to 
correct the data, with pixel edges shown in green. The bottom-right panel shows the scatter in the residuals as a function of binning; the 
red line shows the photon noise limit; the vertical dotted line denotes the timescale of ingress/egress. 
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Fig. 6. — The predicted wavelength-depen dent transit depth of 
WAS P- 12b based on a Solar composition IjBurrows et al.l 120071 . 
I200SI 1. The dotted black line shows a model with a day-like T- 
P profile (large scale height); the solid line shows a model with a 
night-like T-P profile (short scale height). The red points corre- 
spond to a model with ellipsoidal variations; the blue point cor- 
responds to a model without 4.5 fim ellipsoidal variations (see 
first Section 5.3.2). The two black points with error bars o n the 
left show the optical transit depth from Hebb et alj II2009I ) and 
[Macicjewski et al. (2011) (top) and Chan et al. (201]) (bottom); 
we normalize the model transit spectrum to the latter 's observa- 
tion. 
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Fig. 7. — The day-side emergent spectrum of W ASP-12b. The 
color e d line s show various ID atmospheric models HBurrows et al.l 
12007. '2008), while the red points show the measured secondary 
eclipse depths. From l eft to right, the d ata are: z-ban d 
(|l6pcz-Morales eJaDISOTO); J, H, and K^-band llUroll et al.ll20111) . 
IRAC channels 1 & 2 (this study); and IRAC channels 3 & 4 
((Campo et al. 20 1111 . Note that the Kj-band eclipse depth and 
the K-H eclipse color have been confirmed by Zhao ct al. (2011) 
and Crossficld et^] II2012I') . respectively. The Campo ct al. (2011) 
eclipse depths at 3.6 and 4.5 /im are shown in gray. The blue point 
shows the 4.5 fim eclipse depth if ellipsoidal variations arc set to 
zero, the "null hypothesis". In the legend, the first value for 
each model is for the fiducial analysis (including ellipsoidal vari- 
ations), the second value is for the null hypothesis (setting ellip- 
soidal variations to zero). The enhanced-CO model (yellow line) 
offers the best fit, but the Solar composition model (green line) is 
not significantly worse. 

the caveat that these models are in radiative — but not 
chemical — equihbrium, the crtix of fitting WASP-12b's 
unique day-side spectrum is to suppress H2O with re- 
spect to CO. 

Our larger error bars, however, make these composi- 
tion statements marginal: the Solar composition model 



with modest inversion (the green line in Figure [T]) is only 
worse by Ax^ = 8 for 8 degrees of freedom. Furthermore, 
neither the standard composition nor the enhanced CO 
scenario are consistent with the relative transit depths 
at 3.6 and 4.5 /xm (s.f. Figure |6] and previous section). 
It may be possible to reconcile these two measurements 
if the atmospheric composition is grossly different at the 
day-night terminator than near the sub-stellar point. 

5.3. Ellipsoidal Variations 

The best-fit ellipsoidal variations at 3.6 /xm are con- 
sistent with the predicted amplitude of 2 x 10~^, but as 
one can see from the relative uncertainty (or from glanc- 
ing at Figure EJ, they are not robustly detected: the x^, 
residuals and remaining astrophysical parameters do not 
change significantly if ^oiiips is set to zero. 

As shown in Figure 5, however, we clearly detect 
power at 4.5 /im in the second cosine harmonic, cos(2q!), 
consistent with the prediction of ellipso idal variations 
due to the prolate shape of the planet (|Li et al.l 120101 : 
iLeconte et all 120111 : iBudail 1201 If ) . The semi-amplitude 
of the variations, however, is 1-2(2) x 10^'^ using ei- 
ther decorrelation or polynomial IPSV-removal, approx- 
imately 6x the predicted value. 

Since ellipsoidal variations are primarily a geometrical 
effect, it is difficult to understand how the measured am- 
plitude could be be so different at 3.6 and 4.5 /zm. The 
upper layers of the atmosphere should be more distorted; 
the deeper layers more spherical. The relative strengths 
of ellipsoidal variations at the two wavebands imply that 
the 4.5 ^m flux is originating from much higher up in the 
planet's atmosphere than the 3.6 jim flux. The simplest 
way to do this is for the atmosphere to have a greater 
opacity at 4.5 /xm than 3.6 ^m. But the relative transit 
depths indicate exactly the opposite, as described above 
and shown in Figure 

Alternatively, it is possible that detector systematics 
still present after IPSV-removal attenuate the ellipsoidal 
signal at 3.6 /zm or enhance it at 4.5 /xm. The raw pho- 
tometry shown in Figure 1 implies that the 4.5 /im light 
curve is more trustworthy of the two. We therefore begin 
by assuming that the ellipsoidal signal at 4.5 fim is en- 
tirely astrophysical in nature, then consider the opposite 
scenario. 

5.3.1. Interpreting the Ellipsoidal Variations at 4-5 fim 

If we take the 4.5 fim ellipsoidal variations at face 
value, they have surprising implications for the planet's 
shape. The dimensions of a prolate planet may be de- 
scribed by its short and long radii, denoted by Rp and 
^long, respectively. The third dimension (parallel to the 
system's angular momentum vector) is assumed to be 
Rp because we are neglecting rotational effects, which 
tend to produce oblate planets. Note that WASP-12b 
is on a very short-period orbit, so — if tidally locked — 
it has a rotation rate only a factor of two slowe r than 
Jupite r or Saturn. As a result of this rotation, IBudail 
lj201in estimates the planet's polar radius to be 2.5% 
shorter than its lateral equatorial radius^ so strictly 

Note, furthermore, that the sub-stellar and anti-stellar plan- 
etary radii are not equal. Nevertheless, the dominant effect is the 
planet's prolate shape. 
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speaking WASP-12b is a triaxial ellipsoid, but this does 
not affect the analysis below because of the planet's edge- 
on orbit and our relatively short baseline of observa- 
tions is insensitive to t he spin precession of the planet 
(ICarter Winnll2ninbD . 

At conjunction — either transit or eclipse — we are see- 
ing the planet's smallest projected area {TrRp, in the case 
of a perfectly edge-on orbit). The projected area of an el- 
lipsoid on an edge-on orbit varies as fe.g..fVickers 1996i) : 



Ap{a) — ttB1\1 cos^ 



Rr, 



sin^ a. 



(5) 



For Riong/ Rp ~ 1; the changes in projected area follow 
a Ap oc cos(2a) shape and this is in fact how we mod- 
eled them; for more severe elongations, the peaks become 
broader and the troughs narrower, until a limiting case 
of ^„ cx I sinal. 






iime from transit (days) 

Fig. 8. — Our best-fit 4.5 fim phase variations (top panel) can 
be modeled as a geometrical component due to the planet's chang- 
ing projected area (middle panel) and a thermal component due to 
longitudinal variations in the planet's brightness (bottom panel). 
For the red lines all of the 2a power is attributed to the planet's 
ellipsoidal shape; this hypothesis can be ruled out by the transit 
morphology (Figure 9). The blue lines result from assuming that 
the planet is spherical in shape; this hypothesis leads to unphysi- 
cal brightness variations on the planet. The green lines show the 
middle-road: most of the 2a power is due to the planet's shape, 
but the thermal component also contributes a bit. Note that this a 
posteriori analysis has limitations: in our astrophysical model, the 
thermal and ellipsoidal components of phase variations were added, 
while strictly speaking the planet's intensity and cross-sectional 
area should be multiplied. 

If we interpret all of the power in the second cosine 
harmonic as being due to the changing cross-sectional 
area of the planet (the red curves in Figure [8]) , we may 
estimate the planet's aspect ratio as 

day) - 1.8(1). (6) 

The predicted aspect ratio for the planet is R\ong/Rp = 
1.1, while the aspect ratio for a Roche lobe is 3/2 = 1.5. 

It is worth noting that a planet with an aspect ratio 
of 1.8 would have a 13% larger projected area at the 
start and end of transit as compared to transit center, 
resulting in a w-shaped transit, in the absence of stel- 
lar limb-darkening (note that this differs from changes 
in transi t mor phology due to an oblate planet dis- 
cussed in iCarter &: Winn,2010a) . Since the transit depth 



Rlong/ Rp ~ 1 + 2ylcllips(i^*/i^c 



of WASP- 12b is approximately 1%, this shape- induced 
transit effect would come in at the 1.3 x 10"'^ level and 
might be detectable in our current data. In Figure [9] 
we estimate the expected transit morphology by treat- 
ing the planet as a sphere with variable cross-sectional 
area. The red line shows our fiducial model: a spherical 
planet with non-linear limb-darkening of the star. The 
green line shows the effect of the planet's changing cross- 
sectional area (but not its changing shape). The blue line 
shows how limb-darkening partially washes out this sig- 
nal. Figure IHl implies that our data rule out the most 
extreme prolate toy model, but clearly the treatment of 
limb-darkening is important here. 



■ Spherical. w/Lirab- Darkening 

■ Prolate, No Limb-Darkening 

■ Prolate w/Limb- Darkening 




time from transit (days) 

Fig. 9. — WASP-12b transit lightcurve at 4.5 fim; the data have 
been binned for plotting. The red line shows our fiducial model: 
a spherical planet with non-linear limb-darkening of the star. The 
green line shows the effect of the planet's changing cross-sectional 
area (but not its changing shape). The blue line shows how limb- 
darkening partially washes out this signal. For the two prolate 
planet models, we use Rion^/Rp = 1-8, the most extreme scenario 
supported by our phase variations. 

If we instead interpret the phase curve as being caused 
entirely by longitudinal brightness variations on a spheri- 
cal planet (the blue curves in Figure[S]), it can be inverted 
into a longitudinal int ensity map of the planet, following 
iCowan fc Agoll (|2008D . Because of the low-pass filtering 
that occurs in the map — ^ lightcurve convolution, the 
deconvolution will enhance the highest f requency terms 
prese nt in the lightcurve (e.g.. Figure 2 of ICowan fc Agoll 
120091 ). In the current case, the resulting map has two 
prominent temperature peaks: one at the dawn termina- 
tor and one at the dusk terminator. More importantly, 
the only way to simultaneously fit the bright terminators 
and dark night-side is by having negative intensity at the 
anti-stellar point, clearly an unphysical solution. 

Another argument against a spherical planet is that 
while the thermal phase variations of a spherical planet 
will in general contain power in the second harmonic, 
there is no reason for it to all appear in the cosine rather 
than sine term 

When we run fits allowing for the phase 
of the cos(2q!) term to vary, the offset is consistent with 
zero at the 2a level, with the largest offset being less than 
6°. This indicates that the flux is peaking at quadrature, 

'^^ This is in stark contrast to the first harmonic, where one ex- 
pects most of the power to be in the cos a term due to the extreme 
day-night forcing. 
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as expected for ellipsoidal variations. To our knowledge, 
there is no reason to expect such a temperature profile if 
the planet is spherical. 

If the planet is prolate, however, one might expect 
something akin to the gravity darkening/brightening 
seen in some binary stars: the sub-stellar and anti-stellar 
points on the planet are farther from the center of the 
planet than the terminator, leading to lower surface grav- 
ity. The lower surface gravity at the sub-stellar and 
anti-stellar regions might lead to cooler temperatures, 
all things being equal (von Zcipcl 1924). But this would 
only affect the intrinsic component of the planet's power 
budget, expected to be insignificant for hot Jupiters. 

Finally, it is possible that the planet's shape affects the 
circulation of its atmosphere in more subtle ways, leading 
to relatively hot regions at the dawn and dusk termina- 
tors. This brings us to our favored astrophysical interpre- 
tation of the ellipsoidal variations. The planet may have 
an aspect rati oof Rion^/Rp ~ 1-5 (somewhat easing the 
transit morphology constraints described above) , with an 
additional enhancement of the cos(2a) power because of 
a relatively hot day-night terminator (the green curve in 
Figure IH]). Since the atmospheric dynamics on severely 
non-spherical planets has not yet been addressed in the 
literature, it is difficult to say whether this scenario is 
reasonable. 

5.3.2. The Null Hypothesis at 4-5 fJ,m 

If the amplitude of ellipsoidal variations, Adiips, is set 
to zero at 4.5 /xm, the residuals become more correlated 
(as one would guess from Figure 5), increasing the red 
noise and hence uncertainties in the other astrophysical 
parameters. Furthermore, including terms up to sixth- 
order in the polynomial IPSV-rcmoval does not obviate 
the need for the ellipsoidal term, and the same signal 
is detected in the Gaussian decorrelation version of the 
analysis. The is worse when ellipsoidal variations are 
ignored (Ax|; = 0.007 for either the decorrelation or 
polynomial fit). 

However, it is conceivable that the 4.5 /im ellipsoidal 
signal is in fact uncorrected detector systematics. It is 
therefore worth briefly considering the astrophysical im- 
plications of this scenario. 

The most obvious implication of the null hypothesis 
is that the Roche-filling upper-atmosphere of the planet 
need not be optically thick. Since the predicted ampli- 
tude of ellipsoidal variations were only at the 2a level, 
the null hypot h esis is consistent wit h the predic t ions o f 
iLi et all (f2OT0l) . lLeconte"eraII ((2011 and lBudail (fMl . 

Furthermore, the ellipsoidal variations have minima at 
inferior and superior conjunction, so the null hypothe- 
sis causes the transit depth and eclipse depths to sig- 
nificantly increase. Notably, the 4.5 ^m transit depth 
becomes 0.0126(4), comensurate with that at 3.6 ^m. 
The null hypothesis transit depths are consistent with 
a short (night-like) scale height and Solar composition 
(solid black line in Figure HI) . 

The deeper 4.5 ^m eclipse depth, Fday/-F* = 0.0050(4), 
no longer favors the enhan ced CO scenario (yellow line 
in Figure HI) invoked bv Madhusudhan et all (j2011) : the 
fit to the Solar composition model with modest inversion 
(green line in Figure [T]) is only worse by Ax^ ~ 4, for 8 
degrees of freedom. 



5.4. Thermal Phase Variations 

We model thermal phases using only first harmonics 
(cos a and sin a), while our parameterization of ellip- 
soidal variations is a second harmonic (cos 2a). These 
functions are by definition orthogonal, so it is not surpris- 
ing that our conclusions about thermal phase variations 
(^thcim, ttmax) described below are not significantly af- 
fected by the presence or absence of ellipsoidal variations 
discussed above. 

That said, the amplitude and offset we obtain for 
the thermal phase variations depends on which IPSV- 
removal scheme we use. Unfortunately, it is not clear how 
to perform a direct model comparison between decorre- 
lation and polynomial fits using the BIC. The Gaussian 
decorrelation can be thought of as having a large num- 
ber of free parameters and a somewhat smaller number 
of additional constraining equations, but it is not ob- 
vious how to estimate i t s degr ees of freedom. As dis- 
cussed by iBallard et al.l ()2010[ ). the point-by-point de- 
correlation does a great job of removing the short-term 
jitter and long-term detector drift, but may also remove 
any longer-term astrophysical signal. (This did not in- 
terfere with their goal of searching for transits.) Insofar 
as diurnal phase variations are the most gradual astro- 
physical signal in our study, it is not surprising that it is 
the most dependent on IPSV-removal. 

The polynomial fit leads to values at least as good 
— and sometimes significantly better — than the Gaus- 
sian decorrelation. More importantly, the scatter in the 
residuals on the critical 21 minute timescale is lower for 
the polynomial fits, indicating that this method is doing 
a better job of removing correlated noise. 

Our 4.5 ^m pixel sensitivity map obtained from Gaus- 
sian decorrelation (top right panel of Figure 3) shows 
a valley at y « 14.75; this doesn't follow the usual 
pattern of sensitivity decreasing towards pixel edges. 
(Note that we observe the same ripples in sensitivity as 
IB allard et all (HoiO), but those occur on a much smaller 
spatial scale — and exhibit a much smaller amplitude — 
than the y « 14.75 valley.) 

Finally, the middle-left panel of Figure 3 suggests that 
the decorrelation method has over-corrected the 3.6 /xm 
system flux near superior-conjunction: the eclipse bot- 
toms — which should be flat since the planet is hidden 
from view — slope upward towards the central transit. 

We therefore argue that the decorrelation has filtered 
out much of the phase variations along with the system- 
atics and adopt the polynomial-fitted thermal phase pa- 
rameters in what follows Fl 

Following Cowan & Agol (|2011bl ). we estimate the 
hemispheric effective temperatures to be Jday = 
2928(97) K and r„ight = 983(201) K, where the un- 
certainties include an estimate of systematic errors in 
going from brightness temperatures to effective temper- 
aturesl3 This implies a Bond albedo of As — 0.25 
(presumably due to Rayleigh scattering and/or refiec- 

For completeness, we include the thermal phase parameters 
resulting from the decorrelation in Tables 2 and 3. That analysis 
leads to smaller phase amplitudes than the polynomial fit. If taken 
at face value, this implies lower albedo and higher heat transport 
efficiency. 

WASP-12b has a sub-stellar equilibrium temperature of 
To = 3555(132) K, a no-albedo, no-recirculation day-side tem- 
perature of T£=o = (2/3)l/''To = 3213(119) K, and a no-albedo 
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tive clouds) and very low heat recirculation efficiency, 
e < 0.1, at la (see Figure [TUl) PI Note that our ID ra- 
diative transfer models used for interpreting transit and 
eclipse spectra are cloud-free and have an albedo lower 
than this inferred value. 

Space-based optical secondary eclipse depths have been 
measured for a handful of hot Jupiters. If the planet's 
equilibrium temperature is sufficiently low, such mea- 
surements provide an unambiguous es timate of geometri c 
albedo, Ag: 0.04(5) for HP 209458b ('Rowe et al."2008!), 
0.32(3) for Ke pler- 7b jPemorv et a l. 2011), 0.10(2) 
for K epler-17b (IDesert et al.l l2011aD. 0.30(8) for KOI- 
196b (iSanterne et al.M20Tlirand 0.025(7) for TrES-2b 
(jKipping fc Spie gel 2011). Acknowledging that WASP- 
12b is more than 1000 K hotter than any of those planets, 
a Bond albedo of 0.25 is well within the observed range 
for hot Jupiters. 

Assuming gray albedo and a Lambertian scattering 
phase function {Ab = ^Ag-, IHanel et al.l [T992), the 
reflected-light secondary eclipse of WASP-12b should 
have a depth of 2.4 x lO^'', comparable to current ground- 
based precision for this target. In practice, most plan- 
ets exhibit an opposition surge (making them dispropor- 
tionately bright at superior conjunction) and Rayleigh 
scattering, so the actual contrast ratio in blue optical 
wavebands is likely more favorable than this estimate. 




Bond Albedo 

Fig. 10. — 1, 2, and 3a constraints on the Bond albedo and 
recirculation efficiency of WASP-12b from thermal eclipse (blue) 
and phase variatio n (red) observations, using the parametrization 
of Cowan fc Ago] H2011bl ) . The gray-scale shows the confidence 
intervals for the combined constraints. 

Assuming Solar atmospheric composition and hence 
opacity, the 3.6 /xm thermal flux should originate from 
deeper in the atmosphere than any other mid-IR wave- 
band. Insofar as radiative times increase monotonically 
with pressure, we may therefore expect the 3.6 /im phase 
variations to be muted compared to the 4.5 fxni phase 
variations, and the phase offset should be greater at 
the s horter wavelength (e.g.. Figure 9 of iBurrows et al.l 
[20101) . 

fuU-recirculation global te mperature of T^ni = (l/4)^''*To = 
2514(92) K (e.g.. |Cowan~S Agol 2011b). 

If one presumes that advection is the dominant mode of heat 
transport, then e ^ Trad/(Tadv -I- T"rad), where r^^d and r^dv are 
the characteristic radiative and advective timescales at the mid-IR 
photosphere. 



The hot-spot offset is 53(7)° East of the sub-stellar 
point at 3.6 iim. While not as extreme as u- Andromeda b 
(jCrossfield et al.l |2010[ ). it is difficult to reconcile our 
large phase offset and large amplitude at 3.6 /.tm0 It is 
worth noting, however, that there are highly-correlated 
residuals near the purported peak (^ 0.35 day after tran- 
sit) which may be partially responsible for the large off- 
set. 

On the other hand, we find that the 4.5 phase 
amplitude and offset, 16(4)° E, are consistent with a 
iCowan fc: Agoll ()2011aD model with heat transport effi- 
ciency of e = Trad/Tadv « 0.1. To put this in context, 
phase variations and eclip se timing offset o f HD 189733b 
at 8 ^m indicate e « 0.7 (|Agol et al.ll2010[ ). 

5.4.1. Implications of Thermal Phase Variations 

iCowan fc Agoll (|2011bl ) noted that the day-side tem- 
peratures of the hottest short-period giant planets are 
very close to the theoretical upper limit of no albedo 
and no recirculation. This was in contrast to run-of-the- 
miU hot Jupiters (e.g., HD 189733b, HD 209458b), which 
exhibit a variety of albedos/recirculation efficiencies, al- 
beit consistent with generally low albedos {Ab < 0.3). 
In a statistical study of Kepler planetary candidates, 
iCoughlin fc Lopez-Moralei (|2011 [) also found generally 
low albedos for hot Jupiters based on optical secondary 
eclipses. 

The amplitude of the phase variations for WASP- 12b 
depends on the details of the systematics correction, but 
— for reasons stated at the start of Section 5.4 — we favor 
the polynomial fit, which implies a large day-night tem- 
perature contrast, and a non-zero Bond albedo. This 
suggests that the difference between the hottest short- 
period giant planets and other hot Jupiters is not albedo, 
but recirculation efficiency. (Differences in albedo may 
very well explain the differences in day-side effective tem- 
perature amongst the remaining hot Jupiters, however.) 

What could make the hottest of hot Jupiters poor heat 
re-circulators? There are two classes of solutions: de- 
creasing either the planet's characteristic advective fre- 
quency, or radiative timescalePI 

More Magnetic Drag: Assuming these planets have 
magnetic fields, the movement of ionized alkali met- 
als through the field produces drag that is collisionaly 
imparted on the dominant neutral species (presumably 
H and He). Hotter planets should have more ionized 
species, more drag and the refore a ha rder time advect- 
ing heat to their night side (iPerna et~ al. 2010). Because 
of the non-linear dependence of ionization on temper- 
ature, this effect could lead to sudden changes in dy- 
namical regime as one c onsiders increa singly hot planets. 
The scaling relations of iMenoul (|2011l ) indicate that the 
temperature above which magnetic drag severely curtails 
heat transport is inversely related to the planet's mag- 

■^^ It is tempting to attribute the early peak of the 3.6 fim phase 
curve to Doppler beaming, which would produce a peak in flux 
when the star is moving towards us, a quarter period before supe- 
rior conjunction IILoeb fc Gaudi 2003). The expected amplitude of 
this signal, however, is only 7.5 X 10~^ on the Rayleigh-Jeans tail 
of the star. 

One can imagine more exotic means of transporting energy 
(e.g., gravity waves; Watkins & Cho 20l0) but most hydrodynam- 
ical simulations suggest that horizontal energy transport on hot 
Jupiters is primarily a matter of advection. 
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netic field strength. Even the weakest field they consid- 
ered in their study, 3 Gauss, would result in very low 
recirculation efficiency for a planet as hot as WASP-12b. 

Shorter Radiative T imes: Following the argument of 
iCowan &: Agoll (|2011b( ). the radiative rel axation t i me of 
a parcel of gas sc ales as Trad oc (|Iro et all [20051: 
iSeager et alJ 120051 ) . But zonal wind speeds may also 
increase with the amplitude of the diurnal forcing. If 
one assumes that the wind speeds have a fixed Mach 
number, they should scale as Wwind T^/'^, and there- 
fore the advective time should scale as Tadv oc 

r-i/2 (a 

more detailed scaling analysis lead s to the same tem- 
perature dependence; iMenoul [201 1[ ). The stronger de- 
pendence on temperature of radiative time compared to 
advective time implies that — all things being equal — 
hotter planets should be less efficient at balancing their 
day-night temperature contrast. This effect should cause 
the heat transport efficiency to gradually decrease as one 
considers increasingly hot planets. 

Weaker Greenhouse: Atmospheric opacity is typically 
greater at the thermal wavelengths of emergent radiation 
than at the visible wavelengths of incident radiation. For 
the hottest planets, the blackbody peak of thermal emis- 
sion approaches the peak of their host star, so the opac- 
ities of the incoming and outgoing streams converge. In 
that limit, one expects the thermal photosphere and the 
optical deposition depth to be one and the same: thermal 
radiation should escape the atmosphere just as easily as 
the incoming stellar radiation came in. As with the scal- 
ings above, this effect should lead to gradually decreasing 
recirculation efficiency as a function of increasing planet 
temperature. 

More observations of thermal eclipses and phase vari- 
ations for hot Jupiters — especially those near the Tq ~ 
2700 K transition — will be necessary to distinguish 
between the magnetohydrodynamic and the radiative 
timescale arguments. 

6. CONCLUSIONS 

We obtained Warm Spitzer full-orbit phase observa- 
tions of WASP-12b at 3.6 and 4.5 /um, allowing us to 
measure the transit depths, eclipse depths, thermal and 
ellipsoidal phase variations at both wavelengths. We are 
able to push Warm Spitzer photometry to within 10- 
20% of the Poisson limit, but there are two important 
caveats: 

A) Removing intra-pixel sensitivity variations (IPSVs) 
from the data is inherently a model-dependent endeavor. 
This means that we must specify not only an IPSV 
model, but also an astrophysical model before getting 
close to the quoted precision. The simultaneous fit to 
astrophysical and systematic effects makes it difficult to 
produce a "clean" lightcurve independent of astrophysi- 
cal assumptions. For example, we obtain very different 
thermal phase variation parameters depending on how we 
correct for systematics, and it is difficult to distinguish 
between these scenarios based solely on goodness-of-fit. 
Instead, we must resort to a number indirect clues as to 
which IPSV-removal scheme is more trustworthy. 

B) There is still red noise in our residuals, no matter 
how we remove IPSVs. This remaining red noise is the 



dominant source of uncertainty for all of our astrophysi- 
cal parameters. 

We find that WASP- 12b exhibits large-amplitude ther- 
mal phases — indicative of poor day-night heat transport 
and a moderate Bond albedo — but also an unexpectedly 
large phase offset at 3.6 ^m. We do not detect ellipsoidal 
variations at 3.6 /im, while we detect an unexpectedly 
strong signal at 4.5 /im. This leads us to two possible 
hypotheses: 

1) If we take the 4.5 /im ellipsoidal variations at face 
value, we find: deeper transits at 3.6 /im as compared to 
4.5 /im, inconsistent with either Solar or enhanced CO 
models; eclipse depths consistent with previous studies. 
If the 4.5 /im ellipsoidal variations are astrophysical in 
nature, it indicates that the planet is far more distorted 
than predicted, and exhibits a bright terminator. In this 
scenario, the 3.6 /im ellipsoidal variations are attenuated 
due to detector systematics, possibly throwing off the 3.6 
/im transit depth as well. 

2) If instead we presume that the 4.5 /im ellipsoidal 
variations are caused by detector systematics and set 
them to zero — the null hypothesis — we find: transit 
depths consistent with a Solar composition and short at- 
mospheric scale height at the planet's terminator; eclipse 
depths consistent with a Solar composition and a modest 
temperature inversion; ellipsoidal variations in line with 
predictions. 

The null hypothesis is attractive in its simplicity, but 
requires that we were very unlucky; follow-up Warm 
Spitzer observations would have different systematics 
(the PSF would fall on different regions of the pixels) 
and could settle the question of ellipsoidal variations. 
It is likely that near infrared transit spectroscopy could 
break the composition degeneracy, or at least determine 
the atmospheric structure of WASP-12b; if the planet 
has a short scale-height at the terminator it will lend 
credence to the null hypothesis. Further optical transit 
photometry will be useful in pinning down the trans- 
mission spectrum and refining geometrical parameters; 
if a/i?, < 3, then the planet could very well be more dis- 
torted than predicted, making the large ellipsoidal varia- 
tions more plausible. Optical eclipse measurements from 
the ground or from space might confirm the moderate 
albedo of the planet. 

The planet is hypothesized to be losing mass to its 
host star. If this is indeed the case, the presence of an 
accretion disk, accretion stream and impact hot spot may 
necessitate a more holistic model to properly interpret 
observations. 
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